#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <inttypes.h>
#include <regex.h>
#include <math.h>
#include <time.h>

// grater than n
static int n = 6;

// 64 bit ffff
static uint64_t x_ffff = 0xffffffffffffffff; // 0xf*16 = 2^64 - 1

// RNG random number generator
static int p;
static uint64_t s[16];

static inline uint64_t rotl(const uint64_t x, int k) {
	return (x << k) | (x >> (64 - k));}

uint64_t next(void) {
	const int q = p;
	const uint64_t s0 = s[p = (p + 1) & 15];
	uint64_t s15 = s[q];
	const uint64_t result = rotl(s0 * 5, 7) * 9;
	s15 ^= s0;
	s[q] = rotl(s0, 25) ^ s15 ^ (s15 << 27);
	s[p] = rotl(s15, 36);
	return result;}
// end RNG


// rng_init Initialize RNG
static int rng_init(){
    printf("Initializing RNG...\n");
    srand((unsigned)time(NULL)); // seeding rand()
    
    p = (rand()%491)+1; // seeding p, 1<= p <= 491 
    printf("Seeding p with %d\n",p);
    
    for (int i = 0; i < 16; i++) {
        s[i] = rand();
        printf("Seeding s[%d] with 0x%x\n",i,s[i]);} // seeding s[]
    
    int mix = 65537; // 65537, 491
    for (int i = 1; i <= mix; i++) {
        next(); // mixing
        //printf("Initializing XorShift...[%d/%d]\r",i,mix);
        //fflush(stdout);
        }

    printf("\nInitialized RNG successfully.\n");
    return 0;}
// end rng_init


// the_calc Calculation with Theory
static int the_calc_2dice(int num){
    int count_match = 0; 
    //int count_total = 0;
    for (int i = 1; i <= 6; i++){
        for (int j = 1; j <= 6; j++){
            //count_total++;
            if ((i+j>num)&&(j>i)){
                //printf("Got ya!");
                count_match++;}}}
    return count_match;}


// throw_dice Throw a Dice and return a random number from 1 to 6
static int throw_dice(){
    int num;
    //num = (next()/4*3)%6+1;
    num = next()%6+1;
    return num;}


// mc_calc_2dice Calculation with MC 
static uint64_t mc_calc_2dice(int n, uint64_t num){
        int num_dice1, num_dice2;
        uint64_t count = 0;
    for (uint64_t i = 1; i <= num; i++){
        num_dice1 = throw_dice(); //printf("Dice1:%llx\n",num_dice1);
        num_dice2 = throw_dice(); //printf("Dice2:%llx\n",num_dice2);

        if ((num_dice1+num_dice2>n)&&(num_dice1<num_dice2)){
            count++;}}
    return count;}


int main(int argc, char *argv[]){
    printf("Throw 2 Dice Experiment by Z.X.Sun.\n");

    uint64_t throw_num = 0x10000; // 0x10000 = 65536
    if (argc == 2) {
        throw_num = (uint64_t)strtoull(argv[1], NULL, 0);
        printf("Throw 2 Dice %llu times (by user)\n",throw_num);}
    else if (argc == 1){
        printf("Throw 2 Dice %llu times (by default)\n",throw_num);}
    else {
        printf("Usage:    ");
        printf("2dice.run [number]\n");}

    rng_init();

    the_calc_2dice(n);
    int the_count_match = the_calc_2dice(n);
    double the_probabi_2dice = (double)the_count_match/36;
    uint64_t mc_count_match = mc_calc_2dice(n,throw_num);
    double mc_probabi_2dice = (double)mc_count_match / throw_num;

    static double lambda = 1.96; // alpha = 0.05
    static double sigma, error;
    sigma = sqrt( (double)mc_count_match / throw_num
                - pow(((double)mc_count_match / throw_num),2) );
    error = lambda * sigma 
                / sqrt((double)throw_num);

    printf("Theoretical:\n");
    printf("P = %d/36 = %.16f\n",the_count_match,the_probabi_2dice);
    printf("\n");
    printf("Monte Carlo:\n");
    printf("Throw: %llu\n",throw_num);
    printf("Match: %llu\n",mc_count_match);
    printf("P: %.16f\n",mc_probabi_2dice);
    printf("sigma: %.16f\n",sigma);
    printf("error: %.16f\n",error);
    return 0;}